function [selectedFeat selectedFeatIdx] = bciFeatureRanking(runs,nFeat,method,category,validation,doPrint)
% [selectedFeat selectedFeatIdx] =
% bciFeatureRanking(runs,nFeat,method,category,validation,doPrint)
% returns features and corresponding indices of most informative features
% 
% INPUT:
%   runs - runs to involve
%          under consideration: runs can defined as 2 element cell array {bci,Train} to directly
%          apply transformed data  
%   nFeat (optional) - target number of features (if 0<nFeat<1, percentage of 
%                               available number of features is calculated)
%                               default: 0.5
%   method (optional) - string of comparison method:
%             'tval' - t-values (default)
%             'rsqu' - Rsquare values
%             'svm' - weights of a support vector machine
%   category (optional) - which kind of feature
%             'raw' - no reshape of train data, direct use in classifier
%                     space (default)
%             'freq' - select frequencies
%             'band' - select frequencies from normalized predefined bands
%                      (3-10;10-15;15-22;22-30;30-40 Hz) which results in
%                      at least one frequency per band
%             'chan' - select channels
%                      to define the channel index space a prefix can be
%                      prepended:
%                      'good' - good channel space (what classifier sees)
%                      'epoch' - epoch space (the proprocessed channels)
%                      'buffer' - indices in the bvuffer (default)
%                      'actual' - actual channel indices 
%   validation (optional) - validation method to get score over selected
%                           category
%             'max' - maximum value (default)
%             'mean' - mean value
%             'ssq' - sum of squares
%   doPrint (optional) - print results to command window; if category is 'chan' each
%             of the channel index space is printed (default: 0)
% OUTPUT:
%   selectedFeat - best n features 
%   selectedFeatIdx - corresponding indices of category in classification space
%

global GLOBALbci;
global GLOBALtrainDat;

if nargin<3||isempty(method),
    method = 'tval';
end
if nargin < 4 ||isempty(category),
    category = 'raw';
end
if nargin <5 || isempty(validation),
    validation = 'max';
end
if nargin<6,
    doPrint=false;
end

if iscell(runs),
    bci=runs{1};
    Train=runs{2};
    clear runs;
else
    % get the train data
    bci=GLOBALbci;
    % prevent prior feature selection
    bci.param.featureSelection='';
    [bci,Train]=bciGetTrainDat(bci,GLOBALtrainDat{runs});
end

% get the contrasts
if strfind(lower(category),'raw'),
    W = bciMultiCompareFeatures(bci,Train,method,0);
else
    W = bciMultiCompareFeatures(bci,Train,method,1);
end

if strfind(lower(category),'chan'),
    dim = 1; dimT=2;
elseif ~isempty(strfind(lower(category),'freq')) || ~isempty(strfind(lower(category),'band')),
    dim = 2; dimT=1;
elseif strfind(lower(category),'raw'),
    dim = find(size(W{1})>1); 
    dimT=0;
else
    error('Category unknown. Select ''chan'' or ''freq''.');
end

if nargin <2 || isempty(nFeat),
    nFeat = max(length(W),round(0.5*size(W{1},dim)));
end
if nFeat<1,%percentage case 
    nFeat = max(length(W),round(nFeat*size(W{1},dim)));
end

if nFeat>size(W{1},dim),
    error('Number of features available is less than number of features to select');
end

% minimum one feature per class pair
subNfeat = max(1,floor(nFeat/length(W)));
selectedFeatIdx=[];
R2=zeros(size(W{1},dim),length(W));
for k=1:length(W),
    if dimT==0,
        R2(:,k)=abs(W{k});
    elseif strfind(validation,'ssq'),
        % determine sum of squares
        R2(:,k)=sum(W{k}.^2,dimT);
    elseif strfind(validation,'max'),
        % determine maximum
        R2(:,k)=max(abs(W{k}),[],dimT);
    elseif strfind(validation,'mean'),
        % determine mean values
        R2(:,k)=mean(abs(W{k}),dimT);
    else
        error('validation method not defined.');
    end
    % band wise frequency selection
    if strfind(lower(category),'band'),
        % normalise frequencies in predefined bands
        bands=[3 10; 10 15; 15 22; 22 30; 30 40];
        frequencies = bci.param.freq(bci.init.featOfInterest);
        for bandNo = 1: size(bands,1),
            bandIdx=find(frequencies>=bands(bandNo,1)&frequencies<bands(bandNo,2));
            if ~isempty(bandIdx),
                R2(bandIdx,k)=R2(bandIdx,k)*(max(R2(:,k))/max(R2(bandIdx,k)));
            end
        end
    end
    % sort features and get positions
    [R2sorted sortIdx]=sort(R2(:,k),'descend');
    selectedFeatIdx=union(selectedFeatIdx,sortIdx(1:subNfeat));
end
featToAdd = nFeat-length(selectedFeatIdx);
if featToAdd>0,
    unselectedFeatIdx = setdiff(1:size(W{1},dim),selectedFeatIdx);
    [R2sorted sortIdx] = sort(sum(R2(unselectedFeatIdx,:),2),'descend');
    selectedFeatIdx=union(selectedFeatIdx,unselectedFeatIdx(sortIdx(1:featToAdd)));
end

if strfind(lower(category),'chan'),
    goodIdx = selectedFeatIdx;
    epIdx = bciChanMapping(bci,selectedFeatIdx,'good','ep');
    bufIdx = bciChanMapping(bci,selectedFeatIdx,'good','buf'); 
    actIdx = bciChanMapping(bci,selectedFeatIdx,'good','act');
    if strfind(lower(category),'good'),
        selectedFeat = goodIdx;
    elseif strfind(lower(category),'ep'),
        selectedFeat = epIdx;
    elseif strfind(lower(category),'buf'),
        selectedFeat = bufIdx;
    elseif strfind(lower(category),'act'),
        selectedFeat = actIdx;
    else % default
        selectedFeat = bufIdx;
    end
    if doPrint,
%         fprintf('good\tepoch\tbuffer\tactual\n' );
%         for k=1:length(selectedFeatIdx),
%             fprintf('%i\t%i\t%i\t%i\n',goodIdx(k),epIdx(k),bufIdx(k),actIdx(k));
%         end
        disp(selectedFeat);
    end
elseif ~isempty(strfind(lower(category),'freq')) || ~isempty(strfind(lower(category),'band')),
    selectedFeat = bci.param.freq(selectedFeatIdx);
    if doPrint,
        disp(selectedFeat);
    end
elseif strfind(lower(category),'raw'),
    selectedFeat = selectedFeatIdx;
    if doPrint,
        fprintf('%i features selected.\n',length(selectedFeat));
    end
end